library(monocle3)
## Loading required package: Biobase
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## The following object is masked from 'package:Biobase':
## 
##     rowMedians
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## 
## Attaching package: 'monocle3'
## The following objects are masked from 'package:Biobase':
## 
##     exprs, fData, fData<-, pData, pData<-
library(ggplot2)
library(ggpubr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
## 
##     intersect
## The following objects are masked from 'package:IRanges':
## 
##     collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
## 
##     first, intersect, rename, setdiff, setequal, union
## The following object is masked from 'package:matrixStats':
## 
##     count
## The following object is masked from 'package:Biobase':
## 
##     combine
## The following objects are masked from 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.6.2
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(htmlwidgets)
library(NNLM)
library(here)
## here() starts at /data/users/jared/ENS/6mo_LMMP
library(purrr)
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:GenomicRanges':
## 
##     reduce
## The following object is masked from 'package:IRanges':
## 
##     reduce
library(glue)
## 
## Attaching package: 'glue'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## The following object is masked from 'package:SummarizedExperiment':
## 
##     trim
## The following object is masked from 'package:GenomicRanges':
## 
##     trim
## The following objects are masked from 'package:IRanges':
## 
##     collapse, trim
source(here("scripts/accessory_functions/pattern_plotting.R"))
## Loading required package: viridisLite
## 
## Attaching package: 'viridis'
## The following object is masked from 'package:viridisLite':
## 
##     viridis.map
source(here("scripts/accessory_functions/monocle_mods.R"))
source(here("scripts/accessory_functions/geneSetEnrichment.R"))
## 
## clusterProfiler v3.18.1  For help: https://guangchuangyu.github.io/software/clusterProfiler
## 
## If you use clusterProfiler in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.
## 
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:purrr':
## 
##     simplify
## The following object is masked from 'package:IRanges':
## 
##     slice
## The following object is masked from 'package:S4Vectors':
## 
##     rename
## The following object is masked from 'package:stats':
## 
##     filter
## Loading required package: AnnotationDbi
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## DOSE v3.16.0  For help: https://guangchuangyu.github.io/software/DOSE
## 
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:Biobase':
## 
##     combine
## The following object is masked from 'package:BiocGenerics':
## 
##     combine
source(here("scripts/accessory_functions/patternFeatureCorrelationHeatmap.R"))
## ========================================
## circlize version 0.4.12
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================
lmmp_6mo <- readRDS(here("6month_LMMP.rds"))
set.seed(42)

lmmp_6mo <- estimate_size_factors(lmmp_6mo)
print(dim(lmmp_6mo))
## [1] 22794 11123
date <- "10-7"

#Change this to something like cluster, cell type... whatever cell metadata to group by
group_feature <- "clusters"

Generate patterns, or load in existing patterns

do.existing.patterns <- 1
existing_pattern_path <- here("results/NMF/lmmp/old_pattern_run/50dims/")
if(do.existing.patterns){
  date <- ""
  npattern <- 50
  geneWeights.df <- read.csv(paste0(existing_pattern_path, "pattern_gene_weights.csv"), row.names = 1)
  cellWeights.df <- read.csv(paste0(existing_pattern_path, "pattern_cell_weights.csv"), row.names = 1)
} else{

  #use either (size -> log) normalized counts
  #normed_exprs <- normalized_counts(lmmp_6mo, norm_method = "log")
  #or use just log normalized
  normed_exprs <- log10(as.matrix(exprs(lmmp_6mo))+1)
  
  #Choose dimensionality of NMF, how many patterns will be identified
  npattern <- 50
  system.time(decomp <- nnmf(A=as.matrix(normed_exprs),k = npattern,verbose = F, show.warning = T))
  rownames(decomp$H) = paste0("cellPattern",c(1:npattern)) 
  
  ##Gene weights associated with patterns
  geneWeights.df <- as.data.frame(decomp$W)
  colnames(geneWeights.df)<-paste0("cellPattern",c(1:npattern))
  #Merge cell weights with pData for plotting pattern scores on UMAP embedding
  cellWeights.df <- base::as.data.frame(t(decomp$H))
  colnames(cellWeights.df) = paste0("cellPattern",c(1:npattern))
}

#Subset to cell name and feature, group of feature to establish column rows
  cells_by_group <- as.integer(pData(lmmp_6mo)[,group_feature])
  names(cells_by_group) <- colnames(lmmp_6mo)
  #order cells by their cluster, this object is a named integer array
  cells_by_group <- cells_by_group[order(cells_by_group)]
  cells_order <- names(cells_by_group)
  
  
  #Plot heatmap of H (cell weight) matrix
  #Assign a color for each group, these match default ggplot colors
  my_colors <- scales::hue_pal()(length(levels(pData(lmmp_6mo)[,group_feature])))
  names(my_colors) <- unique(cells_by_group)
  ha_cellgroup <- HeatmapAnnotation(cluster = cells_by_group,
                                    col = list(cluster = my_colors))
  
  #We want cells grouped by cluster, as above. We also want to do hierarchical clustering
  #on the pre-defined cluster level. Info must be in the same order 
  cell_mat <- t(cellWeights.df[cells_order,])
  stopifnot(colnames(cell_mat) == cells_order)

Visualize pattern usage by heatmap

plot_cells(lmmp_6mo, color_cells_by = group_feature, group_label_size = 3)
## No trajectory to plot. Has learn_graph() been called yet?
## Warning: The `add` argument of `group_by()` is deprecated as of dplyr 1.0.0.
## Please use the `.add` argument instead.

ComplexHeatmap::Heatmap(cell_mat, name = "cellscore",
                        top_annotation = ha_cellgroup,
                        column_order = cells_order,
                        #cluster_columns = cluster_within_group(cell_mat, groups_to_cluster),
                        show_column_names = F)
## The automatically generated colors map from the 1^st and 99^th of the
## values in the matrix. There are outliers in the matrix whose patterns
## might be hidden by this color mapping. You can manually set the color
## to `col` argument.
## 
## Use `suppressMessages()` to turn off this message.
## `use_raster` is automatically set to TRUE for a matrix with more than
## 2000 columns You can control `use_raster` argument by explicitly
## setting TRUE/FALSE to it.
## 
## Set `ht_opt$message = FALSE` to turn off this message.
## 'magick' package is suggested to install to give better rasterization.
## 
## Set `ht_opt$message = FALSE` to turn off this message.

#row_order = rownames(decomp$H))

Add the pattern weights to pData

#This adds pattern cell weights to pData, make sure they align
#left_join() would be more specific than this cbind()
add.to.cds <- 1
if(add.to.cds){
  stopifnot(rownames(pData(lmmp_6mo)) == rownames(cellWeights.df))
  pData(lmmp_6mo) <- cbind(pData(lmmp_6mo),cellWeights.df)
}

Pattern usage on UMAP

#Color UMAP embedding by cell weights for each pattern.
#Call function and return a list of ggplot objects, and plot to one page
weighted_emb <- lapply(1:npattern,
                       FUN = plotCellPatterns,
                       cds_obj = lmmp_6mo,
                       red_method = "UMAP",
                       do.clip = c(0.02,.98))
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## No trajectory to plot. Has learn_graph() been called yet?
## Cells aren't colored in a way that allows them to be grouped.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
print(weighted_emb[1:npattern])
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

## 
## [[33]]

## 
## [[34]]

## 
## [[35]]

## 
## [[36]]

## 
## [[37]]

## 
## [[38]]

## 
## [[39]]

## 
## [[40]]

## 
## [[41]]

## 
## [[42]]

## 
## [[43]]

## 
## [[44]]

## 
## [[45]]

## 
## [[46]]

## 
## [[47]]

## 
## [[48]]

## 
## [[49]]

## 
## [[50]]

Pattern usage by sample

#Look at pattern usage over sample

condition_patterns<- as.list(1:npattern)
myplots <- lapply(paste0("cellPattern",condition_patterns), 
                  FUN= plotPatternUsageByCondition,
                  cds = lmmp_6mo,
                  bin_by = "sample")
n_by <- 10
seq_by <- seq(0, npattern, n_by)
map(1:(length(seq_by)-1), function(i){
  print(i)
  print(do.call(ggarrange, myplots[(seq_by[i]+1):(seq_by[i+1])]))
})
## [1] 1

## [1] 2

## [1] 3

## [1] 4

## [1] 5

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Save results

if(!do.existing.patterns){
  tmp<-as.data.frame(merge(fData(lmmp_6mo)[,c("gene_id","gene_short_name")],geneWeights.df,by=0))
  dt<- DT::datatable(tmp[,c("gene_id","gene_short_name",
                       unlist(lapply(1:npattern,function(i){paste0("cellPattern",i)}))
  )])
  
  saveWidget(dt, file = paste0("/home/jared/ENS/6mo_LMMP/results/NMF/lmmp/lmmp_6mo_",date,"_pattern_gene_weights.html"))
  
  write.csv(cellWeights.df, paste0("/home/jared/ENS/6mo_LMMP/results/NMF/lmmp/lmmp_6mo_",date,"_pattern_cell_weights.csv"))
  
  write.csv(geneWeights.df, paste0("/home/jared/ENS/6mo_LMMP/results/NMF/lmmp/lmmp_6mo_",date,"_pattern_gene_weights.csv"))
}

Run gene set enrichment analysis

if(!do.existing.patterns){
  #Run GSEA 
  geneWeights.df <- geneWeights.df %>% tibble::rownames_to_column(var = "gene_id")
  
    #putting gene names as rownames "cleaned" them, so genes starting with number or having "-" were changed. 
  geneWeights.df.filt <- geneWeights.df
  geneWeights.df.filt[,"gene_id"] <- rownames(lmmp_6mo)
  
  geneSetEnrichment(gene_weights = geneWeights.df.filt,
                    n_genes = 1000,
                    n_patterns = 50,
                    file_prefix =
                      paste0(here("results/GSEA/lmmp/lmmp_6mo_",date,"_"))
                    )
}

Look at pattern correlations with metadata

#Run pattern correlations
disc_features <- c("cell_type","CCStage", "sample")
cont_features <- c("log10UMI","mito_ratio","tricyclePosition")

pData(lmmp_6mo)$CCStage[is.na(pData(lmmp_6mo)$CCStage)] <- "NA"

heatmaps <- patternFeatureCorrelationHeatmap(lmmp_6mo, cellWeights.df = cellWeights.df, discrete_features = disc_features, continuous_features = cont_features)

hm_edit <- lapply(heatmaps[[1]], function(x){
  x@column_names_param$gp$fontsize <- 8
  x
})

if(!do.existing.patterns){
  
  pdf(here(glue("plots/NMF/lmmp/LMMP_{date}_pattern_feature_correlation.pdf")), width= 12)
  hm_edit[[1]]@row_names_param$gp$fontsize <- 10
  hm_edit
  
  
  
  Heatmap(heatmaps[[2]], column_title = "Correlation of pattern weights with continuous features",
          name = "Pearson",
          column_names_gp = gpar(fontsize = 8),
          column_names_rot = 70,
          row_names_gp = gpar(fontsize = 10))
  
  dev.off()
}

Dont save cds object because pattern names are ambiguous between runs. In the future load the patterns explicitly everytime they are used.